function [idat,gidat] = interp_fast_output(grid,base,data,x)
% [idat,gidat] = interp_fast_output(grid,base,data,x)
%
% Interpolate a time series in an arbitrary point (provided it is
% contained in an element for which the time series is recorded).
%
% grid: the computational grid
% base: finite element basis
% data: the high frequency data structure returned by
%       collect_fast_output
% x: the point where the interpolation is required (column array)
%
% idat:  interpolated data
% gidat: interpolated data gradients
%
% Notice that it is possible to provide different bases for different
% data fields. In this case, the input argument base must be an array
% of bases with one element for each data field.

 d = grid.d;

 % 1) locate the correct element
 toll = 1e-10;
 found = 0;
 ie = 1;
 while( and( not(found) , ie<=length(data.ie) ) )
   e = grid.e(data.ie(ie));
   xi(2:d+1) = e.bi*(x-e.x0); % barycentric coords
   xi(1) = 1 - sum(xi(2:d+1));
   ximin = min(xi);
   if(ximin>=-toll)
     found = 1;
   else
     ie = ie+1;
   end
 end

 if(not(found))
   error(['Can not locate the requested point in any ' ...
          'of the time history elements'])
 end

 % 2) evaluate the basis functions
 if(length(base)==1)
   PHI = zeros(base.pk,1);
   for i=1:base.pk
     PHI(i) = ev_pol(base.p_s{i},xi(2:end)');
   end
   gPHI = zeros(base.me.d,base.pk);
   for i=1:base.pk
     for j=1:base.me.d
       gPHI(j,i) = ev_pol(base.gradp_s{j,i},xi(2:end)');
     end
   end
   % Transform the derivatives in physical space
   gPHI = grid.e(data.ie(ie)).bi'*gPHI;
 else
   for iv=1:length(base)
     PHI{iv} = zeros(base(iv).pk,1);
     for i=1:base(iv).pk
       PHI{iv}(i) = ev_pol(base(iv).p_s{i},xi(2:end)');
     end
     gPHI{iv} = zeros(base(iv).me.d,base(iv).pk);
     for i=1:base(iv).pk
       for j=1:base(iv).me.d
         gPHI{iv}(j,i) = ev_pol(base(iv).gradp_s{j,i},xi(2:end)');
       end
     end
     % Transform the derivatives in physical space
     gPHI{iv} = grid.e(data.ie(ie)).bi'*gPHI{iv};
   end
 end

 % 3) interpolate the time history
 ntimes = length(data.times);
 varnames = fieldnames(data);
 offs = 0;
 for iv=3:size(varnames,1) % skip ie and times
   var_i = getfield(data,varnames{iv}); % i-th variable
   if(ndims(var_i)==3) % add one dimension to scalar fields
     var_i = permute(var_i,[4 1 2 3]);
   endif
   if(length(base)==1)
     P = PHI;
     GP = gPHI;
   else
     P = PHI{iv-2};
     GP = gPHI{iv-2};
   end
   idat (offs+1:offs+size(var_i,1),1:ntimes) = 0;
   gidat(size(GP,1),offs+1:offs+size(var_i,1),1:ntimes) = 0;
   for i=1:size(var_i,1)
     for j=1:length(P)
       idat(offs+i,:) = idat(offs+i,:) + ...
         permute(var_i(i,j,ie,:),[1 4 2 3])*P(j);
     end
   end
   for i=1:size(var_i,1)
     for j=1:size(GP,2)
       for k=1:size(GP,1)
         gidat(k,offs+i,:) = gidat(k,offs+i,:) + ...
           permute(var_i(i,j,ie,:),[1 2 4 3])*GP(k,j);
       end
     end
   end
   offs = offs + size(var_i,1);
 end

%
%   iv = 0;
%
%   nvars = nvars + size(getfield data.)
%
% uup = zeros(d+1,length(data.times));
% for i=1:d
%   for j=1:ubase.pk
%     uup(i,:) = uup(i,:) + ...
%       permute(data.uu(i,j,ie,:),[1 4 2 3])*uPHI(j);
%   end
% end
% for j=1:pbase.pk
%   uup(d+1,:) = uup(d+1,:) + ...
%       permute(data.pp(j,ie,:),[1 3 2])*pPHI(j);
% end

return
